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The Landau level spectrum of graphene superlattices is studied using a tight-binding approach. 
We consider non-interacting particles moving on a hexagonal lattice with an additional one- 
dimensional superlattice made up of periodic square potential barriers, which are oriented along 
the zig-zag or along the arm-chair directions of graphene. In the presence of a perpendicular mag- 
netic field, such systems can be described by a set of one-dimensional tight-binding equations, the 
Harper equations. The qualitative behavior of the energy spectrum with respect to the strength of 
the superlattice potential depends on the relation between the superlattice period and the magnetic 
length. When the potential barriers are oriented along the arm-chair direction of graphene, we find 
for strong magnetic fields that the zeroth Landau level of graphene splits into two well separated 
sublevels, if the width of the barriers is smaller than the magnetic length. In this situation, which 
persists even in the presence of disorder, a plateau with zero Hall conductivity can be observed 
around the Dirac point. This Landau level splitting is a true lattice effect that cannot be obtained 
from the generally used continuum Dirac-fermion model. 
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I. INTRODUCTION 

Recently, the electronic and transport properties of 
graphene superlattices have been the subject of intense 
investigation. Theoretically, it was shown that the pres- 
ence of periodic electrostatioi~— or vector— ~— poten- 
tials, and also of periodic arrays of corrugations 2 ^— 
tailors the graphene properties in a unique way, lead- 
ing to novel features and interesting physics. In one- 
dimensional superlattices, i.e., two-dimensional (2D) su- 
perlattice potentials depending on only one spatial direc- 
tion, the Dirac cones of graphene are distorted, and hence 
the velocity of a particle moving parallel to the potential 
steps is reduced^ Moreover, for certain superlattice pa- 
rameters, this component of the velocity is suppressed, 
and the carriers move only perpendicular to the potential 
steps of the superlattice. Furthermore, for other specific 
superlattice parameters, extra Dirac cones^— and even 
Dirac lines 8 - appear in the energy spectrum of graphene 
besides the usual K or K' Dirac points that exist in the 
continuum model of graphene at the neutrality point in 
the absence of a superlattice. 

Interestingly enough, the emergence of the new Dirac 
points is controlled by the ratio between the potential 
amplitude and the superlattice period, irrespective of the 
superlattice profile, e.g., a cosine^ or a Kromg-Pennej*& 7 - 
type, as long as the period is much larger than graphene's 
lattice constant. The extra Dirac points and their associ- 
ated zero-energy modest drastically affect the transport 
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of the system and also the Landau level 



sequence^ and hence the plateaus in the quantum Hall 
conductivity when a magnetic field is applied. 

The implementation of two-dimensional superlattices 
is another route to modulate the electronic properties of 
graphene. For example, for rectangular superlattices, the 
velocity of carriers is also anisotropic and, depending on 



the Fermi level, the charge carriers are electrons, holes, 
or a mixture of bothi^ Recently, it has been shown that 
for two-dimensional rectangular superlattices the conduc- 
tivity is unchanged from the result of pristine graphene, 
even if the velocity renormalization induced by the su- 
perlattice is quite large Also, new Dirac points with 
and without energy gaps can emerge at high-symmetry 
points in the Brillouin zone in two-dimensional triangu- 
lar superlattices^^ Ab-initio studies of the electronic 
and magnetic properties of graphene-graphane superlat- 
tices have also been reported i 26 ' 27 For instance, it was 
shown that the zig-zag or arm-chair orientation of the 
graphene-graphane interface has a significant impact on 
the electronic properties of the system^ 

Experimentally, there are different possibilities to fab- 
ricate graphene superlattices. For example, it is possi- 
ble to imprint superlattice patterns with periodicity as 
small as 5 nm using electron-beam induced hydrocarbon 
lithography on graphene membranes ^ Graphene grown 
epitaxially on Ru(0001)^ _ — or Ir(lll)^^ surfaces, and 
also on Si C 37 ' 38 shows two-dimensional superlattice pat- 
terns with lattice period of a few nanometers and poten- 
tial strength in the range of few tenths of an electron volt. 
In suspended graphene, the existence of periodic ripples 
has been recently demonstrated^ Another possibility to 
make superlattices with controlled potential amplitude is 
to fabricate periodically patterned gates: p-n and p-n-p 
junctions in graphene^— have already been realized. 

In the present work, we study the evolution of the Lan- 
dau levels of graphene appearing in the presence of a 
one-dimensional superlattice and a strong magnetic field 
applied perpendicular to the graphene plane. For the 
superlattice, we assume a Kronig-Penney type of electro- 
static potential with alternating barriers of +V and — V 
potential strengths and barrier width L a . Here, two cases 
are considered, one with the potential barriers oriented 
along the zig-zag (zz) and the other oriented along the 
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arm-chair (ac) directions of graphene. We use a tight- 
binding approach with nearest-neighbor hopping. The 
magnetic held is incorporated in the hopping integral 
through the Peierls phases as usual. Since the superlat- 
tice depends only on one direction, we use plane waves in 
the other direction. Thus, the two-dimensional problem 
is reduced to a set of one-dimensional equations, known 
as the Harper equations 4£ 

There are two characteristic lengths in the system: one 
is the supcrlatticc barrier width L a and the other is the 
magnetic length I b ■ The behavior of the energy spectrum 
is governed by the relation between L a and Ib- We find 
that if L a is greater than / b , the Landau levels acquire a 
finite broadening (irrespective of disorder) when the su- 
perlattice potential strength increases. In this case, the 
Landau level sequence disappears already for small val- 
ues of V. In the other case, when L a is smaller than Ib, 
the Landau levels picture survives for much higher val- 
ues of the potential strength. Also, two or more Landau 
levels may merge together when V is increased. Conse- 
quently, the plateaus in the Hall conductivity show an 
unconventional sequence that may be controlled by ei- 
ther the external electrostatic potential or the applied 
magnetic field. When the barriers of the supcrlattice po- 
tential arc along the arm-chair direction of graphene, and 
when L a < Ib , we find that the zeroth Landau level splits 
into two sublevels, a novel signature of the interplay be- 
tween the magnetic field and superlattice potential ori- 
entation. We show that the splitting is robust against 
different types of disorder affecting the superlattice po- 
tential. This holds also for a two-dimensional chess-board 
type supcrlatticc. 

The paper is organized as follows. In Section |TT] we 
introduce the theoretical model that we use in order to 
investigate graphene superlattices in the presence of a 
strong magnetic field. In Section IIIII wc present the re- 
sults concerning the energy spectrum of the system for 
various superlattice parameters and magnetic flux den- 
sity strengths. A summary and concluding remarks are 
given in Section HVl 



II. HARPER EQUATION FOR GRAPHENE 
RIBBONS 



The Harper equations of graphene ribbons in a uni- 
form magnetic field reduce the two-dimensional tight- 
binding problem to two one-dimensional equations due 
to the translational invariance along one axis. Deriva 
tions of the Harper equations can be found in Rcf. 
a hexagonal lattice with zig-zag edges and in Rcf. 
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bricklayer lattices with zig-zag and with arm-chair edges. 
Here, we briefly derive Harper equations in the presence 
of a supcrlattice potential. From these, we calculate the 
electronic energy spectra of monolayer graphene ribbons 
with oriented edges in a perpendicular magnetic field and 
a superlattice potential. 

Wc consider a one-band tight-binding model with near- 
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FIG. 1: (a) The profile of the Kronig-Penney superlattice 
potential, (b) Geometry of the zz ribbon with iV Z z zz lines, 
(c) Geometry of the ac ribbon with iV a c dimer lines. Periodic 
boundary conditions are applied in the y direction with length 
L generally large compared to the width W of the sample. 



est neighbor hopping on a hexagonal lattice. The struc- 
ture of graphene ribbons of length L with zig-zag (zz) and 
arm-chair (ac) edges is shown in Fig. [T] The two inter- 
penetrating triangular sublatticcs are denoted by A and 
B. The ribbon width W <C L (x direction) is defined by 
the number N zz of zz lines for the zz ribbon and by the 
number N ar of dimer lines for the ac ribbon 



-N 7 ,,, 



l)a 



=(7V ac -l)^a, 



(1) 
(2) 



where a is the distance between two neighboring carbon 
atoms (a = 1.41 A). 

With ip{ r i) the wavefunction amplitude on site r,, the 
Schrodingcr equation reads 



(3) 



where e is the associated eigenvalue and V the electro- 
static potential (superlattice). The sum runs over the 
nearest neighbors of atom i and Uj is the transfer integral 
between atoms i and j. The magnetic field perpendicu- 
lar to the graphene plane is incorporated in the hopping 
term by means of the Peierls phase 



U 



(4) 



where t is the hopping parameter (t =2.7 eV) and 7^ is 
given by the line integral of the vector potential from site 
Yi to site Tj, and 4>q — h/e is the magnetic flux quantum 
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A-dl. 



(5) 
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We use the Landau gauge A = (0, Bx, 0) and the trans- 
lational invariant direction of each ribbon is taken as the 
y axis. With this particular choice of the gauge, the 
line integral of the vector potential between two nearest 
neighbors i and j becomes 



2tT . X{ -\- Xj 

lij= B{ Vj - yi )^- 



(0) 



The magnetic flux <fi = BS through the area S = 
3a 2 -\/3/2 of a hexagon is taken to be a rational multi- 
ple of the flux quantum <f> = (p/q)(j)Q, hence the magnetic 
flux density B and the magnetic length are set by the 
integers p and q which are chosen to be mutually prime: 



B=^ 



q e 3a 2 V3 



(7) 




The one-dimensional supcrlatticc potential is taken to 
be periodic along the x direction, with periodicity 2L a 
(barrier width L a ) and constant in the y direction 



V( Xi ) =V( Xi + 2L a ). 



(8) 



Figure Q] (a) shows the profile of such a superlattice 
formed by a Kronig-Penney type of electrostatic poten- 
tial. Because the Hamiltonian does not depend on y, we 
use plane waves for the wavefunctions in the y direction 



ip(r) =^(x)e ik * y . 



(9) 



We are ready now to write the Harper equations for the 
zz and the ac ribbon. 



We denote ^f(x^) = ip£, f^fOO = V™ and so the 
Harper equations take the simple form 



B 

m+l 



(12) 



where we defined 



2 cos [k y 
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(13) 



Here, we use = (3m/2— l)a and x^ = (3/2)(m— l)a. 

Equation ()12l) docs not contain any boundary condi- 
tions yet. In the case of zz ribbons with finite N zz , 
Dirichlet boundary conditions are imposed by setting 

V>jv +i = = 0; since there are no atoms for m = 

and m = N zz + 1. Including this fact, Eq. (fi"2"]) simplifies 
for m = 1 and m = N zz to 



(14) 



For periodic boundary conditions (torus), one sets 



ib A/B 



4>\^ B and tpQ 



A/B 



and Eq. (|12p becomes 



(15) 



However, a more careful treatment is needed at the edge 
sites for periodic boundary conditions, due to the require- 
ment for commensurability between the magnetic and the 
lattice lengths: the Peierls phases must be periodic with 
the ribbon width. This implies that 



A. Harper equations for the zz ribbon 

In the following, we label the atoms as shown in 
Fig. [1] (b) and use m to index the zz chains, where m 
takes values between 1 and N zz . Then, the Schrodinger 
equations for the atoms 1A and 2B arc 



e^f(xi) = V(xi)^(xi) + e^e^^Hx*) 



h&ai.B/B 



Vf(<) + ^V>f(x£+i) 



(10) 



^2«) = V(xf n )^(xf n ) + e^e-^VOO 



From Eq. ([5]), the Peierls phases are 

P 1 X m 



7lA,2B = 72B.5A = 27T- 



q 3a 

p 1 X, 
12B,1A = llA,3B = —i-K- 



A i „B 



q 3a 



(ii) 



llAAB = l2B,6A = 0. 



and one finds the condition 

N„ = 



2n Q -, 
P 



(16) 



(17) 



where n can be any nonzero positive integer. This means 
that the parameters for the magnetic flux density (p/q) 
and the ribbon width (given by N zz ) are not independent 
but must be chosen in such a way that Eq. (|T7|) holds. 
In order to obtain the energy spectrum of the system, 
we diagonalize numerically the RHS of Eq. (|12p . i.e., a 
2N ZZ x 2N ZZ matrix with the appropriate boundary con- 
ditions. 



B. Harper equations for the ac ribbon 

To obtain the Harper equations for ac ribbons, we pro- 
ceed in a similar way as in the case of zz edges. Here, m 
indexes now the dimer lines of the ac ribbon, and takes 
values between 1 and -/V ac . The Schrodinger equations for 
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the atoms 1A and 2B of Fig. [TJ(c) are (x^ = xf n = x m ): 

S^(Xrn) = V( Xm )^(x m ) + e**-* ^ V>f (x m ) (18) 

+ e^-'e-^mixm+i) + e^-e-***^ (x m _i) 



Etp2 (Xm.) = V(x m )tp2 ( x m) + e 



H2,i p-ikya I A 



4>1 {x m ) 



with the Peierls phases 

P 2 

llA,2B = ~l2B,\A = 27T- X m (19) 

q 3V3a 

12b.,a = -nw - 2.-^=- 

P 1 X m "4" %m — 1 
72S,6A = -7lA,4B = 271--—= . 

q 3V3a 2 

Again, we denote tpf' B (x m ) = V'm' 5 - Then the Harper 
equations take the simple form 



(20) 



where we defined 

P 2 



a m = cxp |i [k y a + 2tt^ -^—j=x m ] j 
6 m = exp ( - , [k y - + 2tt- — ] } 



(21) 



with x m = (v / 3/2)(m — l)a. 

A/ 13 

For Dirichlet boundary conditions we use tp =0 
and ipff +i — 0) while for periodic boundary conditions 

ipQ B — ip^ B anc ^ ^n^ B +i = V'] 4 I R the case of ac 
ribbons with periodic boundary conditions, commensura- 
bility between the lattice and the magnetic held requires 
that 



A ac = 6n-. 

P 



(22) 



The energy spectrum of an ac ribbon is obtained by diag- 
onalizing the RHS of Eq. ([20l. i.e., a 2A ac x 2A ac matrix. 



III. RESULTS AND DISCUSSIONS 

In this section, we present results obtained by solving 
the Harper equations corresponding to graphcnc super- 
lattices under perpendicular magnetic fields. The super- 
lattice potential is given by a Kronig- Penney function pe- 
riodic along the x direction, with barrier width L a , and 
with alternating ±V barrier heights. The perpendicular 
magnetic field is set by the parameters p and q, according 
to Eq. ((7|) . For simplicity, we take in the following p = 1. 

As a test, we consider first the cases with no supcrlat- 
tice potential or no magnetic field. Then wc show that, 



for V 7^ and fi ^ 0, the energy spectrum and hence 
the system properties strongly depend on the relation be- 
tween the superlattice barrier width L a and the magnetic 
length Z j3 , as well as on the orientation of the superlat- 
tice barriers with respect to the zz or ac directions of 
graphenc. 



A. l/ = 0orB = 

For infinite 2D systems, when the superlattice poten- 
tial is not present (V = 0) and for high magnetic fields, 
the energy eigenvalues are the usual Landau level (LL) 
sequence which, according to the continuum model, oc- 
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18.49 



hv F 

sign(ri) — 1 

IB 



'2nV3p 



n\t, 



(23) 

where n is the LL index, and for the last equality we have 
used Eq. (0 and Hvf = 3/2to. 

The left panel of Fig. [2] (a) shows the energy spectrum 
in the 1 st Brillouin zone of a graphene ribbon in the ab- 
sence of the superlattice potential. The low lying LL en- 
ergies (shown in green, independent of k y ) appear accord- 
ing to Eq. (|23|) . Higher energies (not shown) may deviate 
from the sequence because Eq. ([23| is valid for an infinite 
system. We treat here ribbons with a finite width and are 
only interested in the low-energy domain. For Dirichlet 
boundary conditions, the edge states (k y dependent) are 
shown in red. Between two LLs, the Hall conductivity 
o~ X y is equal to e 2 /hx (number of edge states) and experi- 
ences a jump every time the Fermi energy coincides with 
the LL energy. The corresponding Hall conductivity is 
shown in the right panel of Fig. [2] (a), with the quantized 
values of the integer quantum Hall effect of graphcnei 50 i 51 

In the absence of the magnetic field, and depending 
on the superlattice parameters, the energy spectrum of 
a graphcnc zz superlattice may exhibit additional Dirac 
points for E = 0. For a Kronig-Pcnney superlattice with 
barrier width L a and height V, the number of Dirac 
points increases by two ( x 4 considering valley and pseu- 
dospin degeneracy) whenever the potential amplitude ex- 
ceeds a value of 



U 



Zitta 

i 

2L a 



(24) 



with n a positive integer £ Figure [2] (b) illustrates this 
possibility of changing the number of Dirac points at 
E = by tuning the height of the superlattice poten- 
tial. Shown are the eigenvalues of the Harper equation 
around the Dirac point for a zz ribbon with N zz = 600 
and for periodic boundary conditions in both directions. 
One superlattice barrier contains 60 zz lines, which cor- 
responds to a barrier width of L a = 12.7 nm. One clearly 
sees that when the barrier height V exceeds integer mul- 
tiples of Uq = 0.052 1 (= 0.14cV), then additional zcro- 
encrgy modes appear in the energy spectrum. Our spec- 
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FIG. 2: (a) Left: the energy bandstructure of a graphene zig- 
zag ribbon of width JV ZZ = 600 in the presence of a strong 
magnetic field with p/q = 1/300. The fiat bands correspond 
to the bulk Landau levels. The energies of the k y dependent 
edge states are highlighted in red. Right: The corresponding 
Hall conductivity, (b) The eigenvalues of a graphene Kronig- 
Penney superlattice with barrier width L a = 12.7 nm (60 zz 
lines) for several barrier heights specified in each panel (Uo = 
0.14 eV) in the absence of a magnetic field. 



tra for B = are in agreement with the results obtained 
by means of a continuum Dirac-equation approach for 
graphene superlattices i 5 ' 6 ' 8 ' 10 
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FIG. 3: (a) Energy bands of ac ribbons for V = with 7V ac = 
600 (semiconducting left) and iV a c = 620 (metallic right), (b) 
The energy gap E g for the ac ribbons from (a) in the presence 
of an additional superlattice with L a — 1.22 nm as a function 
of the superlattice potential strength V. 



tor ac ribbon (A^ ac = 600), the energy gap is reduced 
from E g = 0.0162 eV when V = to E g = when 

V = 0.39 eV, and then increases again for higher val- 
ues of the potential strength. For the metallic ac ribbon 
(N^c = 620) the superlattice potential opens a spectral 
gap that reaches a maximum value of E g = 0.022 eV for 

V = 0.45 eV and decreases for higher V. The respective 
minimal and maximal energy gaps E g depend both on 

and L„ . 



C. Superlattice parallel to zig-zag edges 



B. B = and V / ac ribbons 

The electronic properties of graphene nanoribbons 
with ac edges strongly depend on their size. It is known 
that such systems are metallic when (N ac — 2)/3 £ Z 
and semiconductor otherwise^ Figure [3] (a) shows the 
energy bands of two typical ac ribbons, a semiconduct- 
ing one with iV ac = 600 that corresponds to = 73.6 
nm and an energy gap E g = 0.0162 eV, and a metal- 
lic one with = 76.12nm. A superlattice potential 
with L a = 1.22 nm (10 ac dimer lines under each barrier) 
with barriers along the ac edges is imposed on the sys- 
tems and the evolution of the energy gap around E = 
is calculated as a function of the potential strength. The 
results are shown in Fig. [3] (b). For the semiconduc- 



We now consider Kronig-Pcnncy superlattices with 
potential barriers oriented along the zz direction of 
graphene in a strong perpendicular magnetic held. In 
our calculations we use graphene ribbons with N zz = 
12000, which corresponds to a ribbon width of 2555 nm. 
The Kronig-Penney superlattice parameters are chosen in 
such a way that the number of barriers is even. Also, the 
number of magnetic flux quanta per graphene plaquctte 
is fixed at p/q = 1/6000, which gives B = 13.15T for the 
magnetic field and Is = 7.07 nm for the magnetic length. 
These settings allow us to study infinitely long (in the 
y direction) ribbons with Dirichlet or periodic boundary 
conditions in the transverse direction. 

First, we discuss the case when the barrier width is 
larger than the magnetic length. Figure H] shows the re- 
sults for L a = 21.3 nm (100 zz lines per barrier) and 
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FIG. 4: (a) The energy bands of a graphene ribbon with su- 
perlattice with barrier widths of 21.3 nm oriented along the 
zz edges. The applied magnetic field is 13.15 T and the super- 
lattice potential strength is 0.02 eV (left) and 0.05 eV (right). 
The ribbon edge states are highlighted, (b) The energy of 
the LLs as function of the superlattice potential strength, (c) 
The energy-dependent Hall conductivity corresponding to the 
energy bands from (a). 



Ib = 7.07 nm. Typical energy spectra as function of 
the wavevector k y in the first Brillouin zone are given in 
Fig.|4](a) for two different values of the potential strength 
V = 0.02 eV (left) and V = 0.05 eV (right). The differ- 
ence between periodic and Dirichlct boundary conditions 
consists in the appearance of edge states when applying 
the latter. The oscillations seen in the LL energies corre- 
spond to mini-Brillouin zones imposed by the superlat- 
tice. The presence of the superlattice modifies the ener- 
gies of the LLs, which acquire now a finite width depend- 
ing on the value of V. Note that the finite bandwidth of 
the LLs is a consequence of applying a superlattice and 
not because of disorder, which is always present in exper- 
imental situations and induces additional broadening of 
the LLs. The width of the LLs increases with increasing 
V , and for larger values of the potential strength the Lan- 
dau bands merge together in the sense that there remain 
no energy gaps between them. For example, in Fig. [4] (a) 
and for V = 0.05 eV, this is the case for the ±l st and the 
±2 nd LLs. Interestingly enough, the broadening of the 
LLs is not the same for all the LLs. For the parameters 



used in Fig. [4] the width of the ±3 rd LL is smaller than 
that of the other LLs, and this particular LL will merge 
with the others only for higher values of V. The origin of 
this individual broadening is discussed in detail below. 

In Fig. [4] (b) we show the energies of the LLs as a 
function of the superlattice potential strength. The en- 
ergy gaps between the LLs dissappear altogether in this 
case when V > 0.07 eV. Also, for V < 0.07 eV we find 
that the bandwidth of the th LL is equal to the value 
of V. When increasing the superlattice barrier width 
while L a ^> Is, the energy bandwidth of the higher LLs 
also increases with the strength of the potential, and the 
merging of the LLs occurs for even smaller values of V. 

The unusual structure of the energy bands is most di- 
rectly reflected in the plateau sequence of the quantum 
Hall conductivity. Figure @] (c) schematically shows the 
Hall conductivity versus the Fermi energy corresponding 
to the parameters used in panel (a), i.e., L a = 21.3 nm, 
Ib = 7.07 nm. Comparing the cases V = 0.02 eV and 
0.05 eV, one sees that with increasing V the reduction 
of the energy gaps between the LLs leads to a decrease 
of the plateau widths. Moreover, for V = 0.05 eV the 
plateaus at ±(3/2)(4e 2 //i) are not present, and the Hall 
conductivity has an unconventional step size. This is 
because the ±l st and the ±2 nd LLs are now merged to- 
gether, and there is no energy gap between them. Of 
course, the plateaus in the energy- dependent Hall con- 
ductance must not be confused with the experimentally 
observed Hall plateaus, which originate from localization 
due to the intrinsic disorder. 



1. Landau level broadening 

The different broadening of the distinct LLs labeled 
n can be explained using perturbation theory. We con- 
sider the superlattice potential as a perturbation to the 
graphene wavefunctions ip^ belonging to a given energy. 
The integral 

«n(*») = \ J C(*> kvW(x)1%(x, k y )dx (25) 

gives the energy corrections up to first order to one of 
the n th LL energies as a function of k v . The values of 
\n n \ for different LLs provide a quantitative measure of 
the influence of the superlattice on the LL spectrum. The 
energies of the LLs with a larger |k„| are more spread out 
than the energy levels corresponding to a smaller \K n \, as 
illustrated in Fig. [5] Here, we consider a system with 
N zz = 120, p/q = 1/60, V = 0.135 cV and L a = 15 a. 
The energies of the lowest four LLs given in Fig. [5] (a) 
show that the 0th LL (Eq ~ eV) is most broadened by 
the superlattice, and the broadening decreases succes- 
sively for the 1 st LL (E 1 ~ 1.13eV), 2 nd (E 2 ~ 1.56eV) 
and 3 rd (E 3 ~ 1.98 eV). 

The amplitudes of K n reflect directly the broadenings 
of different LLs, as shown in Fig. [5](b). For example, K\ 
has the largest amplitude in comparison with k>i and K3, 
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the latter has the smallest amplitude. Correspondingly, 
in the presence of a superlattice potential, the 1 st LL 
exhibits a larger broadening than the 2 nd , and the 3 rd 
LL shows the smallest broadening. 

To understand the oscillations of the n n integrals with 
k y , we analyze the spatial representation of the wavcfunc- 
tions V'n of graphene in a perpendicular magnetic field. 
Figure [5] (c) shows ^(a;)! 2 for n = 1 (left) and n = 2 
(right) and for two values of k y , which correspond to zero 
(upper part) and maximum (lower part) values of the re- 
spective \K n \ integrals. Note that changing k y only shifts 
the wavcfunctions with respect to the x axis but does not 
alter the structure of l^l 2 . The wavefunctions consist 
of two symmetrical contributions from the A and the B 
sublattices, and the two parts have a reflection symmetry 
axis (dashed line) . The position of the reflection symme- 
try axis of the wavefunctions with respect to the super- 
lattice potential barriers is crucial in determining the am- 
plitude of the respective K n integrals. In the upper parts 
of Fig. M (c), for k y = 2mn/(6aV3),m = 0, 1, 2, . . . the 
symmetry axis of both the 1 st and the 2 nd LLs coincides 
with the superlattice barrier edge, where the potential 
changes sign from +V to — V. In this case, the contribu- 
tion to the K n integral of the wavefunctions situated to 
the left and to the right of the reflection symmetry axis 
cancel each other and k„ takes on a minimal value. On 
the contrary, for k y — (2m + l)n/(6a\/3) the reflection 
symmetry axis coincides with the middle of a superlattice 
potential barrier, the wavefunctions from the left and the 
right sides of the reflection symmetry axis contribute the 
same amount when taking the integral Eq. (|25[) . leading 
to a maximal k„ . For all other fc y , the values of K n fall 
in between these two limiting cases. 

The different amplitudes of n„ for different LLs can 
be also explained from the lower parts of Fig. [5] (c), by 
carefully examining the spatial distribution of the wave- 
function with respect to the barrier width. For the 1 st 
LL, the contributions left and right of the reflection sym- 
metry axis extend mainly over a single superlattice po- 
tential barrier, in this case +V. In the case of the 2 nd 
LL, the wavefunction extends over three barriers. When 
calculating k 2 , one subtracts from the main contribution 
under the +V barrier the part of the wavefunction that 
are extended over the — V barriers. Hence, the amplitude 
of K\ is bigger than the amplitude of K2 and, correspond- 
ingly, the broadening of the 1 st LL is larger than the 
broadening of the 2 nd LL. 

This approach was carried out for several system pa- 
rameters to check its validity. We conclude that such an 
analysis provides an easy way to find out a-priori, start- 
ing only from the unperturbed wavefunctions of graphene 
in a magnetic field, which of the LLs will exhibit a large 
or a small broadening when a superlattice potential is 
switched on. 

When L a rj Ib, the lowest LLs survive for much higher 
values of the superlattice potential strength, as can be 
seen in Fig. [5] (left). Here, L a = 6.4 nm (30 zz lines 
per barrier) and Ib = 7.07 nm. The merging of the th 
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FIG. 5: (Color online) (a) The energies of the four lowest 
LLs of a system with superlattice potential V = 0.135 eV and 
L a = 15 a. The system size and magnetic field are N zz = 120 
and p/q = 1/60, respectively, (b) The k integrals for the 1 st , 
2 nd and 3 rd LLs calculated according to Eq. (|25[) . Here, ki 
(red) has the largest amplitude and K3 (blue) the smallest, 
(c) The wavefunctions of the 1 st (left) and the 2 nd (right) LL 
as a function of the position x for two values of k y where the 
corresponding k„ integrals are zero (upper part) and max- 
imal (lower part). The contributions of the A (blue) and 
the B (green) sublattices to the wavefunction are separately 
shown. Please note that the lines connect only the ampli- 
tudes at the respective lattice sites. The superlattice poten- 
tial is also sketched and the dashed lines are the reflection 
symmetry axes. 



and the ±l st LLs occurs when V takes the value ZUq/2 
(= 0.42 eV for L a = 6.4 nm), which depends only on 
the inverse of the superlattice barrier width according 
to (f2"4l . For V = 3Uo/2 the zero-energy LLs becomes 
three-fold degenerate and, as the Fermi energy is scanned 
from negative to positive energies, the Hall conductivity 
presents a step size of 3(4e 2 //i).— 

The right side of Fig. [6] shows the evolution of the LLs 
as V is increased for a system with L a < Ib- The width 
of the low energy LLs, although finite, is very small, and 
depends weakly on V. A general tendency is the bend- 
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FIG. 6: The evolution of the low-energy LLs as a function of 

V for graphene superlattice along the zz direction and with 
B — 13.15 T (Ib = 7.07 nm). The superlattice barrier width 
is L a = 6.4nm (left) and L a = 2.5 nm (right). 

ing of the LL energies towards E = 0. Again, when 

V = 3U /2 (= 1.05 eV for L a = 2.5 nm, which are 12 zz 
lines per barrier), the th and the ±l 8t LLs merge and 
the zero-energy LL becomes three-fold degenerate. The 
small widths of the LLs in the L a < Ib case can be ex- 
plained using a similar analysis of the n n integrals from 
above. In this case, the wavefunctions of the LLs spread 
over many superlattice barriers. The contributions of the 
wavefunctions under adjacent barriers cancels out when 
calculating the K n integrals, which leads to a very small 
broadening of the LLs. Our results obtained within the 
Harper equation method are consistent with the results 
from a perturbative approach starting from the contin- 
uum model for graphene^ There, it was also found that 
in the case of weak fields (i.e., L a < Ib) the matrix ele- 
ments of the perturbation Hamiltonian do neither depend 
on the center nor on the spread of the LL wavefunctions, 
which leads to flat bands (i.e., small widths of the LLs). 

We have performed several calculations for different 
magnetic held strengths and superlattice parameters, and 
the results shown here are most illustrative. The quali- 
tative behavior of the LL energy spectrum in the three 
regimes, namely L a > Ib, L a ps Ib and L a < Ib is ro- 
bust against changing the system parameters, with the 
prerequisite that both the superlattice barrier width and 
the magnetic length are much larger than the graphene 
lattice constant. 



D. Superlattice parallel to arm-chair edges 

In the following, we discuss the case of graphene su- 
perlattices in a strong magnetic field with the potential 
barriers oriented along the ac direction. When no super- 
lattice is present, the LL energy spectra for ribbons with 
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FIG. 7: The evolution of the low-energy LLs as a function 
of V for graphene superlattices along the ac direction with 
L a > Ib (left) and L a ~ Ib (right). The magnetic flux density 
is B = 13.15 T (Ib = 7.07 nm) and the barrier widths are L a = 
24.6 nm corresponding to 200 dimer lines per barrier (left), 
and L a = 7.3 nm corresponding to 60 dimer lines (right). 



zz or ac edges are identical if periodic boundary condi- 
tions are applied. That is because the LLs are a property 
of the bulk of the graphene ribbon and not of the edges. 
However, when a one-dimensional superlattice potential 
is switched on, then the system can be considered to con- 
sists of many 'sub-ribbons' with internal edges between 
them. Although we still investigate superlattice barrier 
widths much larger than the graphene lattice constant, 
we show below that the orientation of the sub-ribbon 
edges with respect to the zz or ac directions of graphene 
plays an essential role when the magnetic length Ib is 
larger than the sublattice barrier width L a . 

Figure [7] shows the energy of the LLs of ac superlattices 
as a function of the potential strength for two different 
regimes: L a > Ib (left) and L a ps Ib (right). Here, 
p/q = 1/6000 (B = 13.15T and l B = 7.07 nm), L a = 
24.6 nm corresponds to 200 dimer lines per barrier (left) 
and L a = 7.3 nm corresponds to 60 dimer lines per barrier 
(right). The behavior of the LLs is qualitatively similar to 
the case of zz superlattices: for L a > Ib the LLs acquire a 
large broadening and merge together when increasing the 
potential strength, and the LL sequence breaks down for 
even smaller values of V. Also, when L a ps Ib, the LLs 
bend towards zero energy, and their broadening is not so 
strong, so that the LL picture survives for higher values 
of the superlattice potential strength. However, there 
are some significant quantitative differences between the 
ac and the zz cases. For instance, the merging of the th 
and the ±l st LLs does not occur at V = 3Uq/2 any more, 
but at values of V which are more complicated to predict 
from a continuum model, as they do not depend only on 
the inverse of the superlattice barrier width. 

For ac graphene superlattices with strong magnetic 



9 



field, the L a < Ib regime is the most interesting one. In 
this case, the electrons travel, in classical terms, over sev- 
eral potential barriers before closing a cyclotron radius, 
and the ac edge of each 'sub-ribbon' has an unexpected 
influence on the energy spectrum of the LLs. We find a 
splitting of the th LL into two sub-bands when the po- 
tential barrier is increased. Figure [5] (a) illustrates this 
effect for an ac superlatticc with L a = 1.47 nm (12 dimcr 
lines per barrier) and Ib = 7.07 nm. The th LL splits 
as soon as V is turned on, and the energy difference be- 
tween the two subbands increases continuously with V. 
For very large V, a splitting of the higher Landau bands 
was observed too (not shown). Note also that the higher 
LLs bend again towards zero energy, and their broaden- 
ing is very weak. In the presence of Dirichlet boundary 
conditions, no edge states do appear within the energy 
range between the split Landau level. 



E. Origin of Landau level splitting 

A possible explanation for the splitting can be found 
by carefully examining the superlatticc barrier potential 
step. For ac superlattices, the barrier step asymmetri- 
cally divides the graphene hexagons, with 4 atoms on 
one side and two next-neighbor atoms on the other side 
(see the inset of Fig. [H] (a)). In the case of zz super- 
lattices, where the splitting does not occur, the superlat- 
tice barrier symmetrically divides the graphene hexagons, 
with 3 carbon atoms under one barrier and the other 3 
atoms under the next barrier with opposite sign. We 
have considered zz superlattices with an artificially im- 
posed asymmetry, realized by placing an atom from each 
divided hexagon under the adjacent barrier, so that a 4- 
2 asymmetry is created for all hexagons associated with 
a barrier potential step. This configuration is schemat- 
ically shown in the upper inset of Fig. [5] (b). Now, the 
energy spectrum clearly shows the splitting of the th LL 
which is, however, not as strong in magnitude as in the 
case of ac superlattices. Another configuration of an arti- 
ficial superlattice with the barriers oriented along the zz 
direction that divides the step hexagons into 4 atoms on 
one side and a pair of next-neighbor atoms on the other 
side, is shown in the lower inset of Fig. [8] (b). Again, 
the th LL is split into two sub-bands as the strength of 
the potential V is increased, with a splitting magnitude 
larger compared to the previous case. 

These examples point already to the origin of this new 
effect. It is the absence of inversion symmetry due to 
the combined influence of magnetic field and superlatticc 
potential that matters in the ac case. And in the zz sit- 
uations discussed above, the artificially imposed asym- 
metry that is responsible for the splitting destroys the 
inversion symmetry as well. Once more, a closer look 
at the wavefunctions in the absence of the superlattice 
together with the application of a (degenerate) pertur- 
bation theory treatment up to second order in the super- 
lattice potential provides a microscopic explanation for 
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FIG. 8: (a) Energy of the LLs as a function of the superlattice 
potential strength for graphene ac superlattices with L a < Ib- 
The splitting of the th LL is clearly seen. The inset shows the 
barrier edge of the ac-oriented superlattice. (b) The splitting 
of the th LL of zz superlattices with broken symmetry (L| z = 
1.25 nm and Ib = 2.23 nm). The smallest splitting (green, 
labeled 1) corresponds to the system shown in the upper inset, 
and the largest splitting (red, labeled 3) to the lower inset. 
For comparison, the splitting of an ac superlattice with LJ C = 
1.23 nm is also shown (blue data, labeled 2). The insets in 
Fig. [8] (b) schematically show the artificial barrier edges of 
two zz superlattices with broken symmetry that divide the 
graphene hexagons into 2 and 4 carbon atoms. 

the LL splitting. In both the zz and ac situations, the 
Landau levels get broadened in first order perturbation 
theory due to the superlattice potential. A splitting of 
the LL, however, is seen in second order only in the ac 
case because of the lacking inversion symmetry. The lat- 
ter is still present in the zz situation (without an artificial 
symmetry breaking) so that energy levels remain degen- 
erate. Thus, the ac superlattice induced LL splitting is 
the generic case, whereas the particular inversion sym- 
metry present in the perfect zz orientation turns out to 
be only a special situation, probably not met in real sam- 
ples. We conclude that this LL splitting is a true lattice 
effect that was not seen before. The splitting is not to be 
found in the usual continuum model descriptions, since 
there, one normally cannot distinguish between ac or zz 
oriented superlattices. 



F. Splitting in 3 TV dimer samples independent of B 

The splitting for ac superlattices is still present when 
changing the system parameters, i.e., £>, L a , and V, as 
long as we are in the L a < Ib regime. Interestingly 
enough, we find that the splitting of the th LL is maxi- 
mal when L a is chosen such that the number of ac dimcr 
lines under each superlattice barrier is a multiple of 3. 
Moreover, in such cases the energies of the split th LL 
do not change with the magnetic field and depend only 
on V and L a . For other values of L ai when the number 
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of ac dimcrs lines under each barrier is not a multiple of 
3, the splitting is still present, but the energy difference 
between the split sub-bands is one order of magnitude 
smaller than in the previous case and depends on the 
strength of the magnetic field as cx B 2 . 

Experimentally, ideal rectangular superlattices with 
sharp potential jumps can hardly be fabricated. To ver- 
ify that the splitting occurs also when the barrier edges 
are smooth, we have performed calculations considering a 
width w between two adjacent barriers where the super- 
lattice potential changes linearly from +V to — V. We 
find that when increasing w, the energy difference be- 
tween the split subbands of the th LL decreases, but the 
splitting is still present even for w = L a /2. 



G. Influence of disorder 

To make sure that the splitting is not induced by an 
artificial hidden symmetry of the Hamiltonian, we have 
considered different kinds of superlattice disorder in our 
calculations. First, we have studied the case of barrier 
width disorder, by allowing 10% of the barriers to have 
the width L a ± SL a . The results are shown in Fig. [S] 
(a) for a system with Is = 4.08 nm, L a = 1.47 nm, and 

V = 0.25 eV. Increasing the value of SL a from 0.12 nm to 
0.36 nm results in a broadening of the split th LL. When 
the value of SL a is further increased, the broadening over- 
laps the splitting. The next case of disorder considered 
is barrier height disorder. The potential strengths of the 
superlattice barriers are allowed to take random values 
in an interval SV around ±V, as is schematically shown 
in the inset of Fig. |H] (b). Again, we find that this type of 
disorder does not destroy the splitting, and only induces 
an additional broadening of the sublevels, which grows 
with increasing disorder strength. Figure [5] (b) shows 
the energy band in the first Brillouin zone of the split th 
LL for a system with L a = 1.47 nm, Ib = 7.07 nm and 

V = 0.3 eV. Here, the disorder strengths is SV = 0.03, 
0.08 and 0.13 eV, respectively. Shown are the results for 
10 different disorder realizations. 

Finally, we have considered ac superlattices with rough 
edges, as schematically shown in the inset of Fig. (c). 
In this case, the superlattice potential depends on the 
y coordinate, and the Harper equation approach can- 
not be used. Therefore, we have directly diagonalized 
the tight-binding Hamiltonian for graphene with super- 
lattice potential along the ac direction and with rough 
edges. Figure H] (c) shows the energy of the split th LL 
for ac superlattices with rough edges, compared with the 
case of perfect edges. In this case, 1.25% of the atoms 
along one barrier edge are taken at random and forced 
to have a potential strength equal to the one of the next 
barrier. We show the results for 100 disorder realiza- 
tions and conclude that the splitting is still present in 
the spectrum. 

The splitting of the th LL leads to the prediction of a 
& X y = plateau to occur in the quantum Hall conductiv- 
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FIG. 9: (a) The energy spectrum of a graphene ac super- 
lattice with barrier-width disorder. For 10% of the barri- 
ers, the width is L a + 8L a with 8L a = 0.12 nm, 0.24 nm, 
and 0.36 nm. (b) The energy spectrum of a graphene ac su- 
perlattice with barrier-height disorder. Here, the potential 
height varies at random between V — SV and V + SV with 
SV/eV = 0.03,0.08,0.13, and V = 0.3eV. (c) The energy of 
the th LL as function of V for ac superlattices with rough 
edges. 



ity, which may be regarded as a hallmark of ac oriented 
graphene superlattices. Most interesting, the splitting 
occurs only when L a < Ib- Because B ~ lg 2 , this cor- 
responds to magnetic fields that are strong enough to 
have LLs, but low enough to be in the L a < Is regime. 
When increasing B, the magnetic length is reduced and 
the splitting disappears for L a > Ib- Here, we have pre- 
sented the results for L a = 1.5 nm, V ranging between 
and 0.85 eV, and Ib = 7.0 nm which corresponds to 
B = 13.1 T and are realistic parameters. Experimentally 
achievable arc L a between 2 — 10 nm, and Ib = 5 — 20 nm 
(B = 2 — 25 T), and superlattice potential strengths of a 
few tenths of an electronvolt. Thus it should be possible 
to check our findings experimentally 
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IV. SUMMARY 

We have investigated the Landau level structure of sin- 
gle layer graphene in the presence of a one-dimensional 
electrostatic square-potential superlattice. The supcrlat- 
tice modulates the Landau level spectrum in a unique 
way that is most directly reflected in a peculiar band 
broadening, in band bondings, and in an unusual plateau 
sequence of the quantum Hall conductivity. We have 
shown that, depending on the magnitude of the super- 
lattice barrier width L a with respect to the magnetic 
length Ib, the energy band structure of the LLs changes 
dramatically. In general, when L a > Ib, the LLs are 
quickly broadened and merge together for small values 
of the superlattice potential strengths. In the opposite 
regime, when L a < lg, the LLs survive for much higher 
values of V, with a general tendency of the higher-energy 
LLs to bend towards zero energy. 

The orientation of the superlattice barriers with re- 
spect to the zz or ac directions of graphene also plays 
a crucial role when L a < Ib- That is because the su- 
perlattice divides the system into many sub-ribbons of 
width L a that become decoupled when increasing the 
strength of the superlattice potential V. When the mag- 
netic length is larger than the sub-ribbon width, an elec- 
tron travels over many barriers before it closes a cy- 
clotron radius. Therefore, the results differ depending 
on the sub-ribbon edge type the electron encounters at 
each barrier jump. When the sub-ribbons have ac edges, 
we found a novel effect that originates from the interplay 
between graphene's hexagonal lattice and the additional 



superlattice potential barriers and can, therefore, not be 
found in the usual Dirac-fermion continuum model de- 
scription. The new observation is the splitting of the 
zcroth LL which occurs with increasing the superlattice 
potential strength. Alternatively, one can tune the mag- 
netic flux density instead and keep the superlattice po- 
tential strength fixed. This intriguing effect is linked to 
the absence of inversion symmetry in the ac case due 
to the presence of both a superlattice potential and the 
magnetic field. The parameters that we used here are ex- 
perimentally accessible, and the peculiar features of the 
electronic structure may be tested directly in transport 
or optical experiments. 

Finally we mention that according to further calcula- 
tions (results not shown) , a splitting of the zeroth Landau 
level occurs also in the presence of truly two-dimensional 
(chess-board type) superlattices, which remains robust 
even in the presence of additional on-site disorder. This 
observation opens the route to consider the charge den- 
sity fluctuation o 54 ' 55 occurring in real graphene samples 
on a length scale between 20 nm and 30 nm as a natural 
2D-superlattice that replaces the artificial one investi- 
gated in our model calculations. We then suggest that if 
the magnetic field is tuned such that the size of the charge 
'puddles' and the magnetic length had the required ra- 
tio, the resulting gap opening could be responsible for 
an insulating behavior and the diverging resistance at 
the Dirac point with the accompanying a xy = Hall 
plateau as has been observed previously in experiments 
by Checkelsky et al£&2Z 
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